Analytically Solvable Model of Nonlinear Oscillations in a Cold but Viscous and 

Resistive Plasma 
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A method for solving model nonlinear equations describing plasma oscillations in the presence of 
viscosity and resistivity is given. By first going to the Lagrangian variables and then transforming 
the space variable conveniently, the solution in parametric form is obtained. It involves simple 
elementary functions. Our solution includes all known exact solutions for an ideal cold plasma and 
a large class of new ones for a more realistic plasma. A new nonlinear effect is found of splitting of 
the largest density maximum, with a saddle point between the peaks so obtained. The method may 
sometimes be useful where Inverse Scattering fails. 

PACS numbers: 52.30.-q,51.20.+d 
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This paper concerns a fluid description of a plasma and 
should be easily understandable to any Reader familiar 
with fluid dynamics. 

Exact solutions involving realistic situations with non- 
linear effects in plasmas are few and far between. Here 
we present a ID analytically solvable model of nonlin- 
ear oscillations in a two-component viscous and resistive 
plasma. We will obtain generalized plasma waves. 

As is well known to plasma physicists, when viscosity 
and pressure are neglected, a simple Lagrangian coordi- 
nate transformation leads to a known, oscillating, non- 
linear solution [ij. Here, a second coordinate transforma- 
tion defined by the initial density profile will be required. 

We assume that our plasma contains one kind of Z = 1 
ion (protons, deuterons or tritons). They represent a uni- 
form and motionless background for electron fluid oscil- 
lations. The latter are mainly driven by the electric field 
E{x). Electron pressure forces will be neglected. To- 
gether with the assumed infinitely heavy ions this elimi- 
nates ion-acoustic modes. Denoting by no (= const) the 
number density of ions, the basic equations describing 
the electron fluid arc: the continuity equation, momen- 
tum transfer equation and Poisson's equation. 
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where the electron temperature Tg is in energy units, I'e 
is the electron viscosity coefficient and rj is the plasma 
resistivity. We assume that Te is sufficiently small, see 
the following condition ([22|) . so that the first term on the 
right hand side in Eq. ([2]) is negligible as compared to 



the second term. Furthermore, Veixjt) will be modelled 
as 
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This will allow us to solve Eqs. ([T|)-(l3]) analytically. 

Analytical formulas for the electron viscosity coeffi- 
cient i>eix,t) and resistivity r]{x,t) involve several ap- 
proximations, see e.g. As a rule one assumes that 
the plasma is quasineutral (rig ~ Ui) and the distribution 
functions are not far from local Maxwellians. Other ap- 
proximations come from the Chappman-Ensgog method, 
where only first order corrections to local Maxwellians are 
included and terms containing derivatives are neglected. 
Therefore, uncertainty factors of two or three cannot be 
excluded. In a Z = 1 plasma, the relevant formulas given 
in the first reference of take the form: 
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where Xq is a slowly varying Coulomb logarithm (typi- 
cally Aq « 10-20). It can be seen that if rie/no is not 
much different from unity and Te ~ const, then both 
quasineutrality is approximately valid and our modelling 
(|3]) is approximately consistent with Eq. ([5a|l . 

Using Eq. (jl]) and neglecting the electron pressure 
term, we can write Eq. ([2]) in the form 
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The first step to an analytical solution of Eqs. ([T])- 
5|) is to introduce Lagrangian coordinates: xo{x,t), the 
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initial position (at i = 0) of an electron fluid element 
which at time t was at x, and time in the electron fluid 
rest frame, r(= t). The basic transformation equations 
between Eulerian and Lagrangian coordinates are (see i 
for more detail) 

1^ = 1+"^^' + i v^{x„r')dr'. (7) 

The continuity equation ([T]) in the electron fluid rest 
frame is simply 



dx 

rie- — = neo(a;o) = ne{x, t = 0). 
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Therefore, if we introduce an auxiliary variable s related 
to xq by the ID transformation: 
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we obtain 
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Eq. pU]) is equivalent to 
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Using Eqs. ([7]) and pT|) . Eq. ([5]) takes the hnear form 

(12) 
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An important point is that E can also be linearly ex- 
pressed in terms of Ve- Indeed, using 
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which follows from Eqs. ^ and H]), and adding this to 
Ve times Eq. ([3]) we end up with 
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where the right hand side is linear in Vg as promised. 
Eqs. ([H]) and ^ lead to 
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which is a linear partial differential equation for We(s,T) 
with constant coefficients. Solutions of such PDEs are 



any superpositions of the normal modes exp[i(fcs — lot)], 
for which Eq. (|15p leads to the dispersion relation 
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Assuming that k is real and superposing the normal 
modes corresponding to the plus and minus signs in lo 
given by Eq. (I16|) we obtain four real solutions: 



Ve = exp(-Q;fcr)/(tjfer)g(fcs), 



(17) 



where f,g = sin or cos. Our choice will he f = g — sin, 
for which Ve{s,T = 0) = and all three variables x, ^ 
and s will have a common origin (x = ^ = s = 0). 

For each r, Ve given by Eq. (|17p is a periodic func- 
tion of s with wavelength A — 27r/fc. By adding higher 
harmonics, obtained from Eq. (I17p on replacing k —^ nk, 
n = 1,2,..., and multiplying by an amplitude, any so- 
lution periodic in s with wavelength A can be obtained. 
Our choice will be 

oo 

t'e = ^ Ai exp(-a„fcT) sin(cj„fcT) sin(?ifcs). (18) 

n=l 

Our equations and final results take a simple and 
universal form if we introduce dimensionless quantities 
which will be denoted by a bar: 
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Tie = ne/no, E — Eek/{mu!p) 
Thus, Eqs. ®, (Uni) and (Ull) are: 
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Integrating the first two over dx and the last one over df 
and using p^ . we end up with equations which define 
all relevant quantities: x, -E, fie and Ve in terms of s and 
f (= f). Dropping bars for simplicity, the final results for 
the dimensionless quantities become: 
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Ve[s, t) ^ 2_^An exp(-Q!„t) sin(ns), (21d) 



Gn{t) = exp(-a„t) 



sin(u;„i) 
a„ h cos(w„i) 
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where one should replace sin(a-'„t)/ti;„ by t ii ujn — 0. 
One can notice that if n^v + rj > 1, w„(= i|a;„|) becomes 
pure imaginary but Gn(t) remains real (sin(w„f)/cj„ — 
sinh(|w„|i)/|ti;„|, and cos(w„t) = cosh(|ct;„|t). 

Eq. (j21cp is only meaningful if the sum is smaller than 
unity, which imposes a limitation on the An- 

One can eliminate the parameter s from Eqs. (j2ip 
thereby making x and t the independent variables. This 
parameter has to be determined in terms of x and t from 
Eq. (|21ap and used in the remaining Eqs. ((2T|) . While 
numerically this is a simple task, analytical formulas are 
complicated (though for s = or tt, we obtain simply 
s — x). At the same time, the parametric form (j2ip . 
which involves simple elementary functions, can also be 
used to plot E{x,t), ne{x,t) and Ve{x,t) (e.g., by using 
ParametricPlotSD of Mathematica, see Figs.[THl]). 

We can expect our solution (|2ip to be realistic if 
ne{x,t) is not much different from unity. Under this as- 
sumption, the electron pressure term in Eq. ([2|) is neg- 
ligible as compared to the electric field if Tg ^ ^nv^^^ = 
e^noX^/n, equivalent to 

T^ey] < 4.5835 X lO^^noA^ or i/7?< 0.124. (22) 

This indicates that at least one of the coefficients v or rj 
must be much smaller than 0.3. One can prescribe these 
coefficients and one of the parameters Te [eV] , or no and 
determine the second one from Eq. (jl9cp and then find 
A from Eq. (|19dp . 

In Figs. [T] and m we present typical examples for a small 
number of modes included. Various spatially periodic 
structures can be produced. Note the relevance of our 
considerations to a wide range of both laboratory and 
space plasmas. 

The electron density ne{x,t) is an even and periodic 
function of x with wavelength 27r. It can be expanded in 
a Fourier series in x with time dependent Fourier coeffi- 
cients; 




FIG. 1: (color online). Plot of ne{x,t) when only A\ — 0.5 
is nonzero, and v = -q — 10~^. This corresponds to no = 
1.45 X 10^* m"^ and A = 4.32 m, if we assume Te = 10 eV 
and Aq/10 = 2. Damping is practically invisible for the times 
presented. 




3(2;, t) — I + ^ Bm{t) cos{mx). 
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FIG. 2: (color online). Plot of ne{x,t) when only Ai — 0.4 
and A2 — 0.2 are nonzero, u — 0.02 and r; = 10~^. This 
corresponds to no = 1.45 x 10^^ m~'^ and A = 30.6 m, if we 
assume Te = 1 eV and Aq/10 — 2. The viscous damping is 
evident after a single period. 



This along with the first equation in ([^D|) integrated over 
dx leads to the Fourier expansion of E{x,t): 
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And finally, Eqs. ^ and ^ lead to 
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from which Ve{x, t) can be determined. 

The coefficients An and Bm{t) are related by 
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FIG. 3: (color online). Plot of ne{x,t) when only -Bi(O) — 
0.55 is nonzero, = 0.015, rj = lO"'^ and N = 20. This 
corresponds to no = 1.45 x 10^^ and A = 11.2 x 10~^ 

m, if we assume Te — 100 eV and Aq/10 — 2. Note that in 
the presence of even weak viscosity -Bi(O) can exceed 1/2, see 
Eqs. dHJ. 




FIG. 4: (color online). Plot of Tie (a;, t) when only Bi(0) = 0.55 
IS nonzero, u = 0.1 = 10"* and = 20. This corresponds 
to no = 1.45 X 10^'' and A = 1.36 x 10~^ m, if we assume 
Te — 10 eV and Aq/10 = 2. The observed bifurcation of the 
maximum under strong viscosity is a new nonlinear effect. 

These integrals are not expressible in terms of elemen- 
tary functions, but if the sum over m or n is truncated 
at some M or N, one can easily calculate numerically as 
many integrals as needed. If either M — 1 or N = 1, 



An or B„i{t) is expressible in terms of Bessel functions. 
Thus using the identity 3 (p. 185) 

/ cos{ax — z sin x) dx — nJa{z), (27) 

we obtain respectively 

An = -(-l)"+V„(nBi(0)) if B„(0) = for m > 1, 
n 

(28) 

5,„(t) = 2JmimAiGi{t)) if A„ = for n > 1. (29) 

Eqs. (|26p result from standard formulas for Fourier co- 
efficients, for the function ne{s,t — 0)~^ or ne{x,t) as a 
function of x, if the integration variable in these formulas 
is changed from s = x — E to x {n^^ ds ~ dx) or con- 
versely from X — s + E to s [rie dx — ds), see Eqs. (|21ap 
and (fTTj) . 

Eqs. ((29)) and ([23|) - ((25)) present a new solution explic- 
itly given in terms of physical variables x and t. Its plot 
ior ly ^ 1] = 10^^ is shown in Fig. [TJ 

The particular solution (^1]) with An given by (1^ re- 
duces to that of [H if = 7y = though in a different no- 
tation. The known condition i?i(0) < 1/2 follows from 
the identity En=i M^^^) = 2(1^ W (p. 366). For 
Bi{0) = 1/2, becomes zero a.t s — t ~ tt. 

The behavior of this solution for i', 77 > is shown in 
Figs. [3] and [4] In a viscous and resistive plasma -Bi(O) 
can exceed 1/2. Furthermore if v is sufficiently large, 
see Fig. El where = 0.1, a new nonlinear effect can be 
noticed, i.e., the largest density maximum splits in two, 
with a saddle point between the peaks. 
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